Astrophysics and Space
Science, 139 (1987), 1–4
[Erratum: vol. 146,
(No. 1, July 1988), 201]
TRANSFORMATION OF GEOCENTRIC TO GEODETIC
COORDINATES WITHOUT
APPROXIMATIONS
KAZIMIERZ M.
BORKOWSKI
Toruñ
Radio Astronomy Observatory, Nicolaus Copernicus University, Toruñ,
Poland
(Received 27 July, 1987)
| Abstract. An exact and relatively simple analytical transform
of the rectangular coordinates to the geodetic coordinates is presented.
It does not involve any approximation and the accuracy of practical
calculations depends exclusively on the round-off errors. The algorithm is
based on one solution to the quartic equation in tg(45° - y/2), where y is the
parametric (or eccentric) latitude. |
Present-day
astrogeodetic work requires a consistent pair of transforms between geocentric
and geodetic coordinates whose high accuracy holds for any altitude and
latitude. The geocentric coordinates are known to be directly calculable from
geodetic coordinates through the use of simple formulae (e.g. Equations (1)
below). It is not so with the reverse transformation. Until recently there
prevailed a widespread believe that the reverse transformation has not any
elementary solution (e.g. Woolard and Clemence 1966, Heiskanen and Moritz 1967,
Izotov et al. 1974, Lang 1974, Hlibowicki 1981, Taff 1985, Pick 1985, Baranov et
al. 1986, The Astronomical Almanac for the Year 1987 (p. K11)). Thus a
rich variety of different algorithms have been devised to find an approximation
to the geodetic latitude, f, and the vertical hight
(altitude) above a reference ellipsoid, h, given the geocentric coordinates: the
equatorial, r, and the polar, z, component. Paul (1973) seems to have been the
first to demonstrate the existance of exact analytical solutions to the geodetic
coordinates. His algorithm was later considerably improved by Heikkinen (1982)
while Hedgley (1976) presented a solution based on entirely different approach.
The common feature of these three algorithms is solving of a full quartic
(biquadratic) equation in some variable (see also Vaníèek and Krakiwsky 1982).
This is normally tedious work though a final practical procedure may sometimes
be reduced to acceptable number of formulas (e.g. Heikkinen 1982).
Unaware of the existance of the above mentioned publications concerning
the exact analytical solutions I set up to work out one. Now it seems to be
considerably simpler than the three of Paul (1973), Hedgley (1976) and Heikkinen
(1982) both in derivation and in practical implementation. Thus it carries the
potential of becoming a standard in future references. The simplicity of my
solution comes from the fact that of the five coefficients of the quartic
equation only two are different from zero or one (in the absolute value). In the
following we shall preset first a complete derivation and then give rules for
practical use.
The standard formulae relating the geodetic and geocentric coordinates
are:
and
|
z = (Cb2/a2 + h) sin f, |
(1b) |
where
|
C = a2/Ö[(a cos f)2 + (b sin f)2] |
(2) |
is known as the radius of curveture of the reference ellipse at the
direction defined by f, and a and b are the semiaxes of
the ellipse.
Eliminating the h parameter from Equations (1) directly leads to (e.g.
Vaníèek and Krakiwsky, 1982)
|
r – z/tg f = (a2 –
b2)/Ö[a2 + (b tg f)2]. |
(3) |
This is the equation that, under somewhat different representa- tions, served
to derive many of approximations met in practice. It is easy to chack that the
substitution:
|
tg f = a(1 – t2)/(2bt) |
(4) |
transforms Equation (3) into the following quartic equation
|
t4 + 2Et3 + 2Ft – 1 = 0, |
(5) |
where
|
E = [bz – (a2 – b2)]/(ar) |
(6) |
and
|
F = [bz + (a2 – b2)]/(ar). |
(7) |
As is well known, the fourth-degree polynomials can always be solved
using only elementary functions, though these general algorithms often prove
impractical. In our case the simple form of Equation (5) suits well the
procedure known as the Ferrari's solution (e.g. Korn and Korn 1961). It uses one
of the three roots of the qubic resolvent:
where
and
to solve for all four roots of the quartic equation. A real root of Equation
(8) is
|
v = –(Q + ÖD)1/3 – (Q – ÖD)1/3 , |
(11a) |
where
For D < 0 the expression
(11a) is equivalent to
|
v = 2 Ö(–P)
cos{(1/3)arccos[–Q/(–P)3/2]}, |
(11b) |
which avoids the complex
arithmetic. Since the case of D < 0 envelopes points not further than about
45 km from the center of geodetic reference ellipsoids, it is rather unimportant
for practical applications but remains a valuable tool in discussion of the
general properties of the transform.
Now, having calculated v, the four roots of Equation (5) can be found
from
|
t = ±Ö[G2 + (F – vG)/(2G – E)] –
G, |
(13) |
where
|
G = [±Ö(E2 + v) + E]/2. |
(14) |
Of the solutions defined by Equations (13) and (14) all are real for D £ 0, and two are real for D > 0. Any one of the real
solutions can serve to find a distinct geodetic coordinate set through Equations
(4) and (1):
|
f = arctg[a(1 – t2)/(2bt)], |
(15a) |
and e.g.
|
h = (r – at) cosf + (z – b) sinf, |
(15b) |
and all of them satisfy
Equations (1). This completes our derivation, but a few more comments are now in
order.
Among the four solutions there is always one real which falls outside the
range [–90°,+90°], the normal range of latitudes, therefore generally the value
of f of Equation (15a) should be placed in the right
quarter by taking into account relative signs of the numerator and denominator
of the argument of the arctangent function.
The rules can be given to select the single conventional solution, which
corresponds to that value of f which lies in the same
quarter as does the geocentric latitude (arctg z/r), and it can be demonstrated
that there is only one such solution among the four. However, we find it more
practical to restrict the considerations to the almost obvious case of a > b,
and to the northern latitudes (positive or zero z value), which together greatly
simplify the rules: just skip the lower signs of the double signs in Equations
(13) and (14) to take the positive square roots only. On the southern hemisphere
the solutions are of course the same, save the sign reversed, thus the second
restriction does not really limit practical usefulness of our algorithm.
Worth noticing is the fact that the presented algorithm lacks symmetry
about the equatorial plane though, of course, it gives symmetric analytical
results. It means that the numerical results will have different round-off
errors on either side of the equator. This may be indicative of a potential for
improvement of the procedure, however, at this stage it is not clear how it
could be utilized.
The auxiliary parameter t can be shown to be the tangent function of one
half of the complement to 90° of the eccentric latitude (called also the
parametric or reduced latitude), y. An algorithm
exactly analogous to the presented above can be obtained by constructing a
quartic equation in tg(y/2). Such a version exhibit a
singularity at the equator (for z = 0), instead of at the poles as in
our case.
The fact that there are three other solutions besides the desired one
makes practically all iterative procedures sensitive to the starting value for
f. If the geocentric latitude is chosen for the start,
a reasonable choice, there will always be regions on the (r,z) plane (inside the
curve D = 0) such that the iterations will converge to a wrong solution. For
example, the point (r[m],z[m]) = (16000,2000) when referred to the present IAU
ellipsoid may be represented by the following (f[deg],h[m]) coordinates: (69.1546512,–6351904.5),
(–4.3033845,–6362215.0), (–66.8170389,–6355613.9) and (–178.0477051,–6394174.1).
In this example the second solution for geodetic latitude lies much closer to
the geocentric latitude (about +7 degrees) than the first (which is
conventionally the only right).
The described algorithm, as embodied in Equations (6), (7) and (9) through
(15), has been implemented in FORTRAN programming language, and extensively
tested against the forward transform (Equations (1)) for rectangular coordinates
spanning a range of from single meters to more than 500 000 km in both
coordinates. Except for the mentioned singularity at the poles no other
practical limitation was found. The program, in the form of a subroutine named
GEOD,
can be obtained from the author on request.
Acknowledgements
The author wishes to express his
gratitude for financial support from the Nicolaus Copernicus Astronomical Center
of the Polish Academy of Sciences under contract CPBP-01.11.
References
Baranov V.N., Bojko, E.G., Krasnorylov,
I.I., Mashimov, M.M., Plokhov, Yu.V., Urmaev, M.S., and Yashin, S.N.: 1986,
Space Geodesy, Nedra, Moscow, p. 27 (in Russian).
Hedgley, D.R.: 1976, NASA Techn. Rep. R-458.
Heikkinen, M.: 1982, Zeitschrift f. Vermess. 107, 207 (in
German).
Heiskanen, W.A. and Moritz, H.: 1967, Physical Geodesy,
W.H. Freeman and Co., San Francisco, p. 181.
Hlibowicki, R. (ed.): 1981, Higher Geodesy and Geodetic Astronomy,
PWN, Warszawa, p. 273 (in Polish).
Izotov, A.A. (ed.): 1974, Fundamentals of Satellite Geodesy, Nedra,
Moscow (in Russian), p. 34.
Korn, G.A. and Korn, T.M.: 1961, Mathematical Handbook for Scientists and
Engineers, McGraw-Hill, New York, p. 24.
Lang, K.R.: 1974, Astrophysical Formulae, Springer-Verlag, Berlin,
p. 497.
Paul, M.K.: 1973, Bull. Géod., Nouv. Ser., No. 108, 135.
Pick, M.: 1985, Studia Geoph. Geod., 29, 112.
Taff, L.G.: 1985, Celestial Mechanics. A Computational Guide for the
Practitioner, J. Wiley and Sons, New York, Ch. 3.
Urmaev, M.S.: 1981, Orbital Methods of Space Geodesy, Nedra, Moscow,
p. 14 (in Russian).
Vaníèek, P. and Krakiwsky, E.J.: 1982, Geodesy: The Concepts, North
Holland Publ. Co., Amsterdam, p. 324.
Woolard, E.W. and Clemence, G.M.: 1966, Spherical Astronomy, Academic
Press, New York, Ch. 3.